library(MASS)
library(gplots)
library(ggplot2)
library(Hmisc)
library(sciplot)
library(RItools)
library(aod)
library(memisc)
library(lme4)
library(arm)
library(gmodels)
library(dotwhisker)
library(stargazer)
library(psych)


remove(list=ls())
set.seed(44)

setwd("/Users/stephenchaudoin/Dropbox/Kill_Scrapes/replication code for IO/")
d <- read.csv(file="phlsurvey_working_data.csv",head=TRUE,sep=",")

### Note: Figure 4 in the manuscript and Figure I.3 here are the same thing

###
# Balance and Demographics
###
# Table I.1
xBalance(treatment_wod ~ libparty + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data=subset(d, treatment_wod == 1 | treatment_ctrl == 1), report="all")
xBalance(treatment_icccontest ~ libparty + notrightparty +  female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data=subset(d, treatment_icccontest == 1 | treatment_ctrl == 1), report="all")
xBalance(treatment_iccnocontest ~ libparty + notrightparty +  female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data=subset(d, treatment_iccnocontest == 1 | treatment_ctrl == 1), report="all")

# Summary Statistics
describe(d[ , c('treatment_ctrl', 'treatment_wod', 'treatment_iccnocontest', 'treatment_icccontest','libparty','notrightparty','female','cebuano', 'tagalog','catholic','education_postsec','income_lower25th','income_26thto75th','newshours_morethan5','married','metromanilancr')], fast=TRUE)




###
# Analysis
###

## Table I.2; WOD to ICC C; War on Drugs outcome

# 1a: Direct, overall effect of WOD - ICC C on support for WOD
# Full (including speeders)
ols.wod.bin <- lm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.wod.bin)
log.wod.bin <- glm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1), family = "binomial")
summary(log.wod.bin)
ols.wod.num <- lm(supportwod_num ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.wod.num)

ols.wod.bin.ctrll <- lm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.wod.bin.ctrll)
log.wod.bin.ctrll <- glm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1), family = "binomial")
summary(log.wod.bin.ctrll)
ols.wod.num.ctrll <- lm(supportwod_num ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.wod.bin.ctrll)

ols.wod.bin.ctrln <- lm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.wod.bin.ctrln)
log.wod.bin.ctrln <- glm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1), family = "binomial")
summary(log.wod.bin.ctrln)
ols.wod.num.ctrln <- lm(supportwod_num ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.wod.num.ctrln)

# Excluding speeder
ols.wod.bin.ns <- lm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0))
summary(ols.wod.bin.ns)
log.wod.bin.ns <- glm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0), family = "binomial")
summary(log.wod.bin.ns)
ols.wod.num.ns <- lm(supportwod_num ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0))
summary(ols.wod.num.ns)

ols.wod.bin.ctrll.ns <- lm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0))
summary(ols.wod.bin.ctrll.ns)
log.wod.bin.ctrll.ns <- glm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0), family = "binomial")
summary(log.wod.bin.ctrll.ns)
ols.wod.num.ctrll.ns <- lm(supportwod_num ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0))
summary(ols.wod.num.ctrll.ns)

ols.wod.bin.ctrln.ns <- lm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0))
summary(ols.wod.bin.ctrln.ns)
log.wod.bin.ctrln.ns <- glm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0), family = "binomial")
summary(log.wod.bin.ctrln.ns)
ols.wod.num.ctrln.ns <- lm(supportwod_num ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1 & speederorslow == 0))
summary(ols.wod.num.ctrln.ns)


## Figure I.1
# Coef plot for WOD DV, all models, WOD -> ICC C
dwplot(list(ols.wod.bin, ols.wod.num, ols.wod.bin.ctrll, ols.wod.num.ctrll, ols.wod.bin.ctrln, ols.wod.num.ctrln, ols.wod.bin.ns, ols.wod.num.ns, ols.wod.bin.ctrll.ns,  ols.wod.num.ctrll.ns, ols.wod.bin.ctrln.ns, ols.wod.num.ctrln.ns),
	ci = 0.90,
	vline = geom_vline(xintercept = 0), vars_order = c("treatment_icccontest")) +
	    xlab("Coefficient Estimate") + ylab("") + xlim(-0.3,0.02) + 
#	    ggtitle("Effect of ICC Contestation Treatment on Support for WOD") +
    theme(
        plot.title = element_text(face = "bold"),
        legend.position = c(0.007, 0.01),
        legend.justification = c(0, 0),
        legend.background = element_rect(colour = "grey80"),
        legend.title.align = .5
    ) + 
#    ylim(breaks=c(0),labels=c("ICC Cont. Art."))  +
    scale_colour_hue(
        name = "Model",
        labels = c("OLS bin.","OLS num.","bin. OLS w/ctrl Lib.", "num. w/ctrl Lib.","bin. OLS w/ctrl NRP","num. w/ctrl NRP","bin. OLS no speed","num. no speed","bin. OLS w/ctrl Lib., no speed","num. w/ctrl Lib., no speed","bin. OLS w/ctrl NRP, no speed","num. w/ctrl NRP, no speed")
    )

# Used for labelling stars in the legend
stargazer(ols.wod.bin, log.wod.bin, ols.wod.num, ols.wod.bin.ctrll, log.wod.bin.ctrll, ols.wod.num.ctrll, ols.wod.bin.ctrln, log.wod.bin.ctrln, ols.wod.num.ctrln, ols.wod.bin.ns, log.wod.bin.ns, ols.wod.num.ns, ols.wod.bin.ctrll.ns, log.wod.bin.ctrll.ns,  ols.wod.num.ctrll.ns, ols.wod.bin.ctrln.ns, log.wod.bin.ctrln.ns, ols.wod.num.ctrln.ns,
		type = "latex",
          title = "Effect of ICC Cont. on Support for WOD",
          label="tab:wodtoicccwod",
          omit.stat=c("f", "ser","adj.rsq"))



## Figure I.2

# 3. Going from ICC NC to ICC C

#a. Direct effect on WOD
# Full (including speeders)
ols.wod.bin <- lm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.bin)
log.wod.bin <- glm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1), family = "binomial")
summary(log.wod.bin)
ols.wod.num <- lm(supportwod_num ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.num)

ols.wod.bin.ctrll <- lm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.bin.ctrll)
log.wod.bin.ctrll <- glm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1), family = "binomial")
summary(log.wod.bin.ctrll)
ols.wod.num.ctrll <- lm(supportwod_num ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.bin.ctrll)

ols.wod.bin.ctrln <- lm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.bin.ctrln)
log.wod.bin.ctrln <- glm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1), family = "binomial")
summary(log.wod.bin.ctrln)
ols.wod.num.ctrln <- lm(supportwod_num ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.num.ctrln)

# Excluding speeder
ols.wod.bin.ns <- lm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.bin.ns)
log.wod.bin.ns <- glm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0), family = "binomial")
summary(log.wod.bin.ns)
ols.wod.num.ns <- lm(supportwod_num ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.num.ns)

ols.wod.bin.ctrll.ns <- lm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.bin.ctrll.ns)
log.wod.bin.ctrll.ns <- glm(supportwod_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0), family = "binomial")
summary(log.wod.bin.ctrll.ns)
ols.wod.num.ctrll.ns <- lm(supportwod_num ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.num.ctrll.ns)

ols.wod.bin.ctrln.ns <- lm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.bin.ctrln.ns)
log.wod.bin.ctrln.ns <- glm(supportwod_bin ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0), family = "binomial")
summary(log.wod.bin.ctrln.ns)
ols.wod.num.ctrln.ns <- lm(supportwod_num ~ treatment_icccontest + notrightparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.num.ctrln.ns)


### Coef plot for WOD DV, all models, ICC NC -> ICC C
dwplot(list(ols.wod.bin, ols.wod.num, ols.wod.bin.ctrll, ols.wod.num.ctrll, ols.wod.bin.ctrln, ols.wod.num.ctrln, ols.wod.bin.ns, ols.wod.num.ns, ols.wod.bin.ctrll.ns,  ols.wod.num.ctrll.ns, ols.wod.bin.ctrln.ns, ols.wod.num.ctrln.ns),
       ci = 0.90,
       vline = geom_vline(xintercept = 0), vars_order = c("treatment_icccontest")) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.25,0.15) + 
  #	    ggtitle("Effect of ICC Cont. on Support for WOD (Base: ICC No Contest)") +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = c(-0.001, 0.01),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title.align = .5
  ) + 
  #    ylim(breaks=c("treatment_icccontest"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS bin.","OLS num.","bin. OLS w/ctrl Lib.", "num. w/ctrl Lib.","bin. OLS w/ctrl NRP","num. w/ctrl NRP","bin. OLS no speed","num. no speed","bin. OLS w/ctrl Lib., no speed","num. w/ctrl Lib., no speed","bin. OLS w/ctrl NRP, no speed","num. w/ctrl NRP, no speed")
  )

# Regression table, used for labelling stars in the figure legend
stargazer(ols.wod.bin, log.wod.bin, ols.wod.num, ols.wod.bin.ctrll, log.wod.bin.ctrll, ols.wod.num.ctrll, ols.wod.bin.ctrln, log.wod.bin.ctrln, ols.wod.num.ctrln, ols.wod.bin.ns, log.wod.bin.ns, ols.wod.num.ns, ols.wod.bin.ctrll.ns, log.wod.bin.ctrll.ns,  ols.wod.num.ctrll.ns, ols.wod.bin.ctrln.ns, log.wod.bin.ctrln.ns, ols.wod.num.ctrln.ns,
          type = "latex",
          title = "Effect of ICC Cont. on Support for WOD (Base: ICC No Contest)",
          label="tab:iccnctoicccwod",
          omit.stat=c("f", "ser","adj.rsq"))




## Figure 4 in main manuscript / Figure I.3 in appendix
fig.reg <- lm(supportwod_bin ~ treatment_wod * libparty + treatment_icccontest * libparty + treatment_iccnocontest * libparty, data = d)

# Calculating predicted values from model
fig.dataset.temp1 <- data.frame(treatment_wod=c(0,1,0,1), treatment_icccontest=c(0,0,0,0), treatment_iccnocontest=c(0,0,0,0), libparty=c(0,0,1,1))
fig.dataset.temp2 <- data.frame(treatment_wod=c(0,0,0,0), treatment_icccontest=c(0,1,0,1), treatment_iccnocontest=c(0,0,0,0), libparty=c(0,0,1,1))
fig.dataset.temp3 <- data.frame(treatment_wod=c(0,0,0,0), treatment_icccontest=c(0,0,0,0), treatment_iccnocontest=c(0,1,0,1), libparty=c(0,0,1,1))

fig.dataset <- rbind(fig.dataset.temp1, fig.dataset.temp2, fig.dataset.temp3)
fig.fitted <- predict(fig.reg, fig.dataset, type=c("response"), se.fit=T)

# Creating the plot
coef <- c(fig.fitted$fit[1], fig.fitted$fit[7], fig.fitted$fit[2], fig.fitted$fit[4], fig.fitted$fit[10], fig.fitted$fit[12], fig.fitted$fit[6], fig.fitted$fit[8])
#1 control nonlib # 7 control lib ### 2 wod nonlib # 4 wod lib ### 6 icccontest nonlib # 8 icccontest lib ### 10 iccnocontest nonlib # 12 iccnocontest lib

se <- c(fig.fitted$se.fit[1], fig.fitted$se.fit[7], fig.fitted$se.fit[2],fig.fitted$se.fit[4], fig.fitted$se.fit[10],fig.fitted$se.fit[12], fig.fitted$se.fit[6], fig.fitted$se.fit[8])
lb <- coef-1.96*se
ub <- coef+1.96*se
pid <- c(1.00,1.00,0,0,-1.00,-1.00,-2.00,-2.00)
var <- c("Control","Control","WOD article","WOD article","ICC NC article","ICC NC article","ICC C article","ICC C article")
party <- c("NL","L","NL","L","NL","L","NL","L")

results <- data.frame(coef,se,lb,ub,pid,var,party)

col <- c("GREY","BLACK","GREY","BLACK","GREY","BLACK","GREY","BLACK")
pd <- position_dodge(width=0.01)

supportforwod <- ggplot(results, aes(coef,pid,color=col,alpha=.5)) +
  geom_point(aes(shape=party),size=4) + 
  geom_segment(aes(x = lb, y = pid, xend = ub, yend = pid, linetype=c("1","1","1","1","1","1","1","1")),position=pd, lwd=1.5) +
  scale_color_manual(values=c("grey40","black")) + 
  scale_shape_manual(values=c(19,18)) +
  scale_x_continuous("Predicted support for War on Drugs", limits=c(0.38,0.93),breaks=c(0.40,0.50,0.60,0.70,0.80,0.90),labels=c("0.40","0.50","0.60","0.70","0.80","0.90")) + 
  scale_y_continuous("",limits=c(-2.25,1.35),breaks=c(-2,-1,0,1),labels=c("ICC Contest Art.","ICC No Contest Art.","WOD Article","Control")) +
  theme_bw(base_size = 12) +
  theme(panel.grid.minor = element_blank()) +
  theme(panel.grid.major = element_blank()) +
  theme(axis.ticks.x=element_blank()) +
  theme(axis.ticks.y=element_blank()) +
  theme(legend.position = "none")  +
  theme(axis.text.x=element_text(size=11)) +
  theme(axis.text.y=element_text(size=12)) +
  theme(legend.title=element_blank()) +
  theme(legend.key = element_blank()) +
  annotate("text", x = 0.58, y = 0.15, label = "Liberal Party") +
  annotate("text", x = 0.78, y = 0.15, label = "Non-Liberal Party") +
  annotate("text", x = 0.88, y = 0.85, label = "0.88") +
  annotate("text", x = 0.79, y = 0.85, label = "0.79") +
  annotate("text", x = 0.58, y = -0.15, label = "0.58") +
  annotate("text", x = 0.78, y = -0.15, label = "0.78") +
  annotate("text", x = 0.77, y = -1.15, label = "0.77") +
  annotate("text", x = 0.45, y = -1.15, label = "0.45") +
  annotate("text", x = 0.75, y = -2.15, label = "0.75") +
  annotate("text", x = 0.53, y = -2.15, label = "0.53")

supportforwod

# Regressions with different comparisons
fig.reg2 <- lm(supportwod_bin ~ treatment_icccontest * libparty, data = subset(d, treatment_iccnocontest==1 | treatment_icccontest==1))
fig.reg2 <- lm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_wod==1 | treatment_icccontest==1))
fig.reg2 <- lm(supportwod_bin ~ treatment_icccontest * libparty, data = subset(d, treatment_wod==1 | treatment_icccontest==1))

# Regression for the claim "The net result is that support levels are similar to those observed among respondents who read a ``regular'' WOD article."
fig.reg2 <- lm(supportwod_bin ~ treatment_icccontest, data = subset(d, treatment_wod==1 | treatment_icccontest==1))




## Table I.3, regressions for support for WOD, by treatment, party

ols.wod.binl <- lm(supportwod_bin ~ treatment_icccontest * libparty, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.binl)
ols.wod.numl <- lm(supportwod_num ~ treatment_icccontest * libparty, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.numl)
ols.wod.bin.ctrll <- lm(supportwod_bin ~ treatment_icccontest * libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.bin.ctrll)
ols.wod.num.ctrll <- lm(supportwod_num ~ treatment_icccontest * libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1))
summary(ols.wod.num.ctrll)

ols.wod.binl.ns <- lm(supportwod_bin ~ treatment_icccontest * libparty, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.binl.ns)
ols.wod.numl.ns <- lm(supportwod_num ~ treatment_icccontest * libparty, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.numl.ns)
ols.wod.bin.ctrll.ns <- lm(supportwod_bin ~ treatment_icccontest * libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.bin.ctrll.ns)
ols.wod.num.ctrll.ns <- lm(supportwod_num ~ treatment_icccontest * libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_wod != 1 & speederorslow == 0))
summary(ols.wod.num.ctrll.ns)


## Table I.4
# Effect of ICC C treatment on support for ICC
ols.icc.bin <- lm(iccsupport_bin ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.icc.bin)
ols.icc.num <- lm(iccsupport_num ~ treatment_icccontest, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.icc.num)

ols.icc.bin.ctrl <- lm(iccsupport_bin ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.icc.bin.ctrl)
ols.icc.num.ctrl <- lm(iccsupport_num ~ treatment_icccontest + libparty + female + cebuano + tagalog + catholic + education_postsec + income_lower25th + income_26thto75th + newshours_morethan5 + married + metromanilancr, data = subset(d, treatment_ctrl != 1 & treatment_iccnocontest != 1))
summary(ols.icc.num.ctrl)

# Short table, (null) effect of WOD -> ICC C on support for the ICC, in aggregate
stargazer(ols.icc.bin, ols.icc.num, ols.icc.bin.ctrl, ols.icc.num.ctrl,
          type = "latex",
          title = "Effect of ICC Cont. on Support for ICC",
          label="tab:wodtoicccicc",
          omit.stat=c("f", "ser","adj.rsq"))

## Figure I.4
# Predicted support for ICC, party interactions

fig.reg <- lm(iccsupport_bin ~ treatment_wod * libparty + treatment_icccontest * libparty + treatment_iccnocontest * libparty, data = d)

# Calculating predicted values from model
fig.dataset.temp1 <- data.frame(treatment_wod=c(0,1,0,1), treatment_icccontest=c(0,0,0,0), treatment_iccnocontest=c(0,0,0,0), libparty=c(0,0,1,1))
fig.dataset.temp2 <- data.frame(treatment_wod=c(0,0,0,0), treatment_icccontest=c(0,1,0,1), treatment_iccnocontest=c(0,0,0,0), libparty=c(0,0,1,1))
fig.dataset.temp3 <- data.frame(treatment_wod=c(0,0,0,0), treatment_icccontest=c(0,0,0,0), treatment_iccnocontest=c(0,1,0,1), libparty=c(0,0,1,1))

fig.dataset <- rbind(fig.dataset.temp1, fig.dataset.temp2, fig.dataset.temp3)
fig.fitted <- predict(fig.reg, fig.dataset, type=c("response"), se.fit=T)


# Creating the plot
coef <- c(fig.fitted$fit[1], fig.fitted$fit[7], fig.fitted$fit[2], fig.fitted$fit[4], fig.fitted$fit[10], fig.fitted$fit[12], fig.fitted$fit[6], fig.fitted$fit[8])
#1 control nonlib # 7 control lib ### 2 wod nonlib # 4 wod lib ### 6 icccontest nonlib # 8 icccontest lib ### 10 iccnocontest nonlib # 12 iccnocontest lib

se <- c(fig.fitted$se.fit[1], fig.fitted$se.fit[7], fig.fitted$se.fit[2],fig.fitted$se.fit[4], fig.fitted$se.fit[10],fig.fitted$se.fit[12], fig.fitted$se.fit[6], fig.fitted$se.fit[8])
lb <- coef-1.96*se
ub <- coef+1.96*se
pid <- c(1.02,.90,0,0,-1.00,-1.00,-2.00,-2.00)
var <- c("Control","Control","WOD article","WOD article","ICC NC article","ICC NC article","ICC C article","ICC C article")
party <- c("NL","L","NL","L","NL","L","NL","L")

results <- data.frame(coef,se,lb,ub,pid,var,party)

col <- c("GREY","BLACK","GREY","BLACK","GREY","BLACK","GREY","BLACK")
pd <- position_dodge(width=0.01)

supportforicc <- ggplot(results, aes(coef,pid,color=col,alpha=.5)) +
  geom_point(aes(shape=party),size=4) + 
  geom_segment(aes(x = lb, y = pid, xend = ub, yend = pid, linetype=c("1","1","1","1","1","1","1","1")),position=pd, lwd=1.5) +
  scale_color_manual(values=c("grey40","black")) + 
  scale_shape_manual(values=c(19,18)) +
  scale_x_continuous("Predicted support for ICC", limits=c(0.50,1.00),breaks=c(0.60,0.70,0.80,0.90),labels=c("0.60","0.70","0.80","0.90")) + 
  scale_y_continuous("",limits=c(-2.25,1.35),breaks=c(-2,-1,0,1),labels=c("ICC Contest Art.","ICC No Contest Art.","WOD Article","Control")) +
  theme_bw(base_size = 12) +
  theme(panel.grid.minor = element_blank()) +
  theme(panel.grid.major = element_blank()) +
  theme(axis.ticks.x=element_blank()) +
  theme(axis.ticks.y=element_blank()) +
  theme(legend.position = "none")  +
  theme(axis.text.x=element_text(size=11)) +
  theme(axis.text.y=element_text(size=12)) +
  theme(legend.title=element_blank()) +
  theme(legend.key = element_blank()) +
  annotate("text", x = 0.68, y = 0.15, label = "Non-Liberal Party") +
  annotate("text", x = 0.85, y = 0.15, label = "Liberal Party")

supportforicc




